#Matching Analysis
load("Amnesty Provision Main Data.Rdata")
#Labels
labs<-c("Incompatibility",
        "Civil Society Strength",
        "Electoral Democracy",
        "Political Violence",
        "Rule of Law")

#Estimate Genetic Matching######
match<-na.omit(data[,c("Amn", "bin", "Inc", "v2x_cspart", "v2x_polyarchy", "v2x_clphy", "v2x_rule")])
treat_formula<-as.formula("bin~Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule")
set.seed(1015)
m.out<-matchit(treat_formula, data=match,
               method="genetic",
               replace=F,
               pop.size=1000) #eventual pop.size=1000 if genetic

#Table A5 - Balance Statistics
balance<-data.frame(summary(m.out)[5])
balance<-balance[2:6,"reduction.Mean.Diff."]
balance<-cbind.data.frame(labs, balance)
balance[,2]<-round(balance[,2],2)
colnames(balance)<-c("Variable",
                     "/% Improvement in Difference in Means")
balance
balance<-print(xtable(balance, digits=3, align=c("l", "l", "c")),
               include.rownames = F)
write(balance, file="Table 5.tex")

#Figure 8
colnames(m.out$X)<-labs
plot(m.out, type="QQ", interactive=F)

pdf("Fig8.pdf")
print(plot(m.out, type="QQ", interactive=F))
dev.off()


#Stochastic Simulation for Figure 5#####
set.seed(1015)
mdat<-match.data(m.out)
zelig<-zelig(Amn~bin+Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule, data=mdat,
             model="relogit")

#Probability of Amnesty w/o Substitutes
sim.no<-setx(zelig, bin=0,
             Inc=1,
             v2x_cspart=mean(mdat$v2x_cspart),
             v2x_polyarchy=mean(mdat$v2x_polyarchy),
             v2x_clphy=mean(mdat$v2x_clphy),
             v2x_rule=mean(mdat$v2x_rule))

#And with Substitutes
sim.bin<-setx(zelig, bin=1,
              Inc=1,
              v2x_cspart=mean(mdat$v2x_cspart),
              v2x_polyarchy=mean(mdat$v2x_polyarchy),
              v2x_clphy=mean(mdat$v2x_clphy),
              v2x_rule=mean(mdat$v2x_rule))

#Get the first difference
sim.out<-sim(zelig, x=sim.no, x1=sim.bin)

#Gather the plot info
ev.no<- sim.out$get_qi(qi ="ev", xvalue = "x")
ev.bin<- sim.out$get_qi(qi ="ev", xvalue = "x1")

plot_dat<-cbind.data.frame(ev.no, ev.bin, ev.bin-ev.no)
names(plot_dat)<-c("1) No Substitutes", "2) Substitutes", "3) Difference")

#Get observed CI
gglines<-cbind.data.frame(apply(plot_dat,2,sort)[25,],
                          apply(plot_dat,2,sort)[975,])
gglines$variable<-cbind(rownames(gglines))

gglines<-reshape2::melt(gglines, id="variable")
names(gglines)[2]<-"open"

plot_dat<-reshape2::melt(plot_dat)

p<-ggplot(plot_dat, aes(value, color=variable, fill=variable))+
  geom_density(alpha=.75, color=NA) + 
  labs(color="Amnesty Provision", fill="Amnesty Provision")+
  theme(legend.position="top")+
  geom_vline(data = gglines,aes(xintercept = value, color=variable), size=.25) +
  geom_vline(aes(xintercept=0), color="black", size=.25) +
  ylab("Density")+
  xlab("Simulated Probability of Peace Agreement Amnesty")+
  facet_grid(variable~.)+
  theme_bw()+
  scale_fill_manual(values=viridis_pal()(4)[1:3])+
  scale_color_manual(values=viridis_pal()(4)[1:3])+
  theme(legend.position="none")

p
pdf("Fig5.pdf",  width=6, height=4)
print(p)
dev.off()


#Bounds Analysis#####
set.seed(1015)
Tr<-match$bin
X<-match[,!colnames(match)%in%c("bin", "Amn")]
Y<-match$Amn

gen1 <- GenMatch(Tr=Tr, X=X, pop.size=1000,replace=FALSE)


mgen1 <- Match(Y=Y, Tr=Tr, X=X, Weight.matrix=gen1, replace=FALSE)

summary(mgen1)

MatchBalance(bin~Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule,
             data=match, match.out=mgen1, nboots=1000)

#Rosenbaum Bounds#####
binarysens(mgen1, Gamma=2, GammaInc=.1)

#Treatment model for comparison
#Want to scale the inputs for interpretability
match[,3:7]<-scale(match[,3:7], center=T)
treat<-brglm(bin~Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule, data=match)

#Average and maximum impact on treatment
median(exp(coef(treat)))
max(exp(coef(treat)))

#Topline number from paper- ratio of hypothetical to meidan
#Observed effect
1.5/median(exp(coef(treat)))


#Estimate Full Propensity Score Matching #####
match<-na.omit(data[,c("Amn", "bin", "Inc", "v2x_cspart", "v2x_polyarchy", "v2x_clphy", "v2x_rule")])
treat_formula<-as.formula("bin~Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+v2x_rule")
set.seed(1015)
m.out.full<-matchit(treat_formula, data=match,
                    method="full",
                    replace=F, distance = "logit")

#Table A6 - Balance Statistics
balance<-data.frame(summary(m.out.full)[5])
balance<-balance[2:6,"reduction.Mean.Diff."]
balance<-cbind.data.frame(labs, balance)
balance[,2]<-round(balance[,2],2)
colnames(balance)<-c("Variable",
                     "/% Improvement in Difference in Means")
balance
balance<-print(xtable(balance, digits=3, align=c("l", "l", "c")),
               include.rownames = F)
write(balance, file="Table 6.tex")

#Figure 9
colnames(m.out.full$X)<-labs
plot(m.out.full, type="QQ", interactive=F)

pdf("Fig9.pdf")
print(plot(m.out.full, type="QQ", interactive=F))
dev.off()


mdat<-match.data(m.out.full)

#Stochastic Simulation for Figure A10
set.seed(1015)
zelig<-zelig(Amn~bin+Inc+v2x_cspart+v2x_polyarchy+v2x_clphy+ v2x_rule, data=mdat,
             model="relogit")
summary(zelig)

#Probability of Amnesty w/o Substitutes
sim.no<-setx(zelig, bin=0,
             Inc=1,
             v2x_cspart=mean(mdat$v2x_cspart),
             v2x_polyarchy=mean(mdat$v2x_polyarchy),
             v2x_clphy=mean(mdat$v2x_clphy),
             v2x_rule=mean(mdat$v2x_rule))

#And with Substitutes
sim.bin<-setx(zelig, bin=1,
              Inc=1,
              v2x_cspart=mean(mdat$v2x_cspart),
              v2x_polyarchy=mean(mdat$v2x_polyarchy),
              v2x_clphy=mean(mdat$v2x_clphy),
              v2x_rule=mean(mdat$v2x_rule))

#Estimate first difference

sim.out<-sim(zelig, x=sim.no, x1=sim.bin)

ev.no<- sim.out$get_qi(qi ="ev", xvalue = "x")
ev.bin<- sim.out$get_qi(qi ="ev", xvalue = "x1")

plot_dat<-cbind.data.frame(ev.no, ev.bin, ev.bin-ev.no)
names(plot_dat)<-c("1) No Substitutes", "2) Substitutes", "3) Difference")

gglines<-cbind.data.frame(apply(plot_dat,2,sort)[25,],
                          apply(plot_dat,2,sort)[975,])
gglines$variable<-cbind(rownames(gglines))

gglines<-reshape2::melt(gglines, id="variable")
names(gglines)[2]<-"open"

plot_dat<-reshape2::melt(plot_dat)

p<-ggplot(plot_dat, aes(value, color=variable, fill=variable))+
  geom_density(alpha=.75, color=NA) + 
  labs(color="Amnesty Provision", fill="Amnesty Provision")+
  theme(legend.position="top")+
  geom_vline(data = gglines,aes(xintercept = value, color=variable), size=.25) +
  geom_vline(aes(xintercept=0), color="black", size=.25) +
  ylab("Density")+
  xlab("Simulated Probability of Peace Agreement Amnesty")+
  facet_grid(variable~.)+
  theme_bw()+
  scale_fill_manual(values=viridis_pal()(4)[1:3])+
  scale_color_manual(values=viridis_pal()(4)[1:3])+
  theme(legend.position="none")

p
pdf("Fig10.pdf",  width=6, height=4)
print(p)
dev.off()

rm(list=setdiff(ls(), c("cluster_bootstrap", "int_plot")))
